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Abstract 

We consider a transformed Ornstein-Uhlenbeck process model that can be 
a good candidate for modelling real-life processes characterized by a combina- 
tion of time-reverting behaviour with heavy distribution tails. We begin with 
presenting the results of an exploratory statistical analysis of the log prices of 
a major Australian public company demonstrating several key features typical 
of such time series. Motivated by these findings, we suggest a simple trans- 
formed Ornstein-Uhlenbeck process model and analyze its properties showing 
that the model is capable of replicating our empirical findings. We also discuss 
three different estimators for the drift coefficient in the underlying (unobserv- 
able) Ornstein-Uhlenbeck process which is the key descriptor of dependence in 
the process. 

Keywords: Ornstein-Uhlenbeck process, heavy tails, regular variation, rank cor- 
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1 Introduction 

There exist a large number of mathematical models designed to reproduce the dynamics 
of financial data time series and aimed at capturing the key features of their behaviour, 
one of the most important of them being heavy distribution tails. A recent monographic 
reference presenting an up-to-date overview of the field is [2], two further monographic 
references devoted specifically to various aspects of using heavy-tailed distributions 
being [231 H] - Despite the abundance of (sometimes rather sophisticated) such models, 
one cannot claim that there is single widely accepted satisfactory model (or even class 
of models) for the behaviour of the log prices of traded financial assets. Among the 
main reasons for that is the ever-changing economical environment that must have 
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substantial effect on the dynamics of the financial time series and hence keeps making 
the already verified and fitted models obsolete. Furthermore, different uses (which 
may include options pricing and/or forecasting, at different time scales) of the models 
stipulate different requirements on their structure and properties, which means there 
can hardly be a single answer to all the requests. Hence the continuing interest in 
considering alternative mathematical models that may be successfully used to either 
get insights in the nature of the log price dynamics or even solve some practical applied 
problems. 

The present note may be considered as a follow-up to a relatively recent paper [IT] 
which analysed a class of rather simple diffusion models for the log prices. The class 
consists of transformed Ornstein-Uhlenbeck (OU) processes of the form Y t = h(X t ), 
where h is a smooth enough strictly increasing function, 

dX t = -aX t dt + rcLWt, t > 0, (1) 

Wt is the standard Wiener process, and a, r > are constants. Clearly, both transition 
probabilities and stationary distribution are readily available for the process Y t , and it 
was shown in [11] that this class of models (i) allows a closed-form expression for the 
likelihood function of discrete time observations, (ii) allows the possibility of heavy- 
tailed observations and (iii) allows the analysis of the tails of the increments. The 
paper also discussed fitting the model to the log share prices, via approximating by it 
the hyperbolic diffusion process considered in [4j , where the latter model was found to 
possesses a number of properties encountered in empirical studies of stock prices and 
was "rather successfully fitted to two different (Danish) stock price data sets" . 

However, despite one of the listed in [TT] motivations for considering the class being 
the need to be able to model heavy tailed data, the paper only goes as far as to show 
(Proposition 2) how to choose h to obtain tail behavior of the form 

P(y t >x) = x c exp(-ax 5 + 6)(1 + O^ 1 )), x — y oo, 

with 5 G [1, 2], a > and b, c £ R, while Theorem 3 of the paper gives conditions on 
h under which the transition density will have exponential decay at infinity. This kind 
of tails can hardly be classified as being heavy. In fact, as it is well-known and will be 
confirmed in our Section [2] below, one often deals with a power function tail decay for 
real-life data sets. 

The main objectives of the present paper include eliciting conditions on h under 
which the distribution tails (for the stationary law of h(X t ), the distributions of the 
increments h(X t ) — h(X s ) and the transition probabilities) for the process Y t = h(X t ) 
will have regularly varying behaviour at infinity thus agreeing with empirical observa- 
tions, and showing that the dependence structure of such a model is also consistent 
with a typical Australian stock log price behaviour. We used BHP Billiton Ltd prices 
from 2003-2010 for illustration purposes, but obtained similar results when doing our 
calculations for several other major companies, both from mining and other sectors 
of the economy as well, that are listed on the Australian Securities Exchange (ASX), 
including the ANZ Banking Group Ltd and TOLL Holdings Ltd. 

Further arguments in favour of the proposed simple model include the fact that there 
is good Gaussian fit for the empirical copulas of the process' values, and the observation 
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that the parameter a characterising the dependence structure of the process can be 
easily estimated using rank correlations techniques, and that these estimates computed 
from the time series increments for different time lags agree with each other. 

One of the important aspects of mathematical modelling of time series is the choice 
of the time scale which the modelling will be aimed at. Firstly, the dynamics at different 
scales can be different, as the key driving factors at these scales may be different — and 
changing because of changing environment, as it may be happening with the small 
scale picture over the last years, due to the wide spread of fast algorithmic trading. 
Secondly, capturing the detail of the small-scale behaviour of prices can be of little 
help if one is interested in medium or long-term modelling. In our study, we are mostly 
interested in medium-term modelling, for daily data with time horizons of the order 
of 10 — 50 days. The empirical data we used seem to confirm the appropriateness of 
the suggested model: even though some of its parameters may be changing with time, 
it appears that its structure remains consistent with observations. 

For further reading and references concerning using OU processes for modelling of 
financial time series, the interested reader is referred to [19]. 

The paper is organised as follows. In Section [2] we present the results of an ex- 
ploratory analysis of the BHP Billiton Ltd log prices time series, to elicit the key 
features that successful candidate models (for time horizons around 10 days) must 
possess (tail behaviour, dependence structure). Section [3] presents a theoretical anal- 
ysis of the model Y t = h(X t ), {X t } following (II]), demonstrating conditions on h that 
lead to regularly varying distribution tail behaviour and also showing that the process' 
dependence structure (in terms of rank correlations) is consistent with the real life data 
properties discussed in Section [2j Section [5] contains a number of comments (mostly 
concerned with the estimation of a) on fitting the suggested model to data. Section [5] 
presents the proofs of the theorems formulated in the paper. 

2 Distributional properties of log returns 

A lot of work has been done on studying the distributional properties of the stock (log) 
price time series, and we refer the interested reader to [21 HI [231 E] and the numerous 
references therein for illuminating discussions of the past findings. However, for the 
purposes of this research, we decided to do new statistical analysis of a single stock 
price time series, with a view of looking at the aspects of the empirical data set that 
can help us clarify the suitability of the transformed OU process model. In this section, 
we use the time series of daily (closing) prices St for BHP Billiton Ltd stock as quoted 
by the ASX, t = 0,1,..., N = 2061, the day t = being 1 January 2003, t = N 
corresponding to 17 December 2010, to elicit the typical distributional properties of 
the log prices Z t := \nS t . As we pointed out earlier, we did the same calculations for 
the stock prices of several other major Australian companies from different sectors of 
economy (e.g. ANZ Banking Group Ltd and TOLL Holdings Ltd), which yielded quite 
similar results, so it appears that the data used in this section may be deemed to be 
typical. 

All the data analysis and graphics were done in Matlab. 
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2.1 The tail behaviour 



The usual starting point in data analysis of this kind is to make normal Q-Q-plots. 
Since we are interested in a dynamic model for the data, we do that for the increments 

A s Z t := Z t - Z t -&, t = 5, 5 + 1,..., N, 

for different values of 5. Clearly, A$Z t is the log return of the stock over a period of 
time of length 5 ending at time t. A typical result is shown in Fig. [T] displaying the plot 
for 5 = 20: there is a very good straight line fit in the central part (representing more 
than 95% of all data, so that the bulk "middle part" of the distribution looks pretty 
much normal), and then in the end regions there are clear deviations from the straight 
line indicating the presence of heavy tails. Different time lag values in the range up to 
5 = 100 return rather similar pictures. 
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Figure 1: The normal Q-Q plot for the 20 days log-returns {A 20 Zt}. 

To elicit the character of the tail behaviour of the empirical distributions, one can 
use the log-log plots of their tails. Figure [2] displays the natural logarithms of the 
empirical distributions' left and right tails for the log returns {AsZ t } against the logs 
of the data values, for four different time lag values: 5 = 1,5,10 and 20. The key 
common features of the plots are: a nice parabolic shape in the "middle region" (which 
indicates a good normal fit in that region), and then a very good straight line fit for 
the "remote regions" (comprising about 5% of all data) in all the cases, indicating the 
power decay of the distributions. An interesting observation is that, for each of the right 
tails, the slope of the fitting straight line (and hence the index of the approximating 
power function) is one and the same, for all time lag values. 

The numerical (absolute) values j3±(5) of the slope coefficients of the fitted straight 
lines for different values of 5 are given in Table [I] (the subscripts +/— correspond to the 
right/left tails). The table also shows the values of the Hill estimators 0±,{ P }(S) for the 
exponents of the power functions specifying the tail decay rate [10]. Recall that, under 
the assumption that the right tail of the theoretical distribution of the data V\, . . . , V n 
with order statistics V\, n > ••• > V n ,n has a regularly varying form x~@ + L(x), L(x) 
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Figure 2: The log-log plots of the empirical distributions' tails of the BHP stock price log 
returns for time periods of 5 = 1, 5, 10, 20 (days). The left (right) pane presents the plots for 
the left (right, resp.) tails. Straight lines are fitted in each case using the mean squares (the 
total number of observations from these regions comprise about 5% of all data). 

being a slowly varying function as x — > oo, the Hill estimator /3 +) { p } of {3 + is computed 
according to 

k— 1 

1 1 v — -\ Vj f, k 
- = -, t> ln-^, p=100-. 

Note that the estimator was shown to be consistent (as n — > oo, p — > 0) for observations 
forming a strongly mixing sequence [25] . 

Observe that the estimates obtained by both methods are generally in good agree- 
ment, especially for the right tail which appears to be "heavier". Note that for that 
tail the values of the estimates for the regular variation index remain roughly the same 
for all 5 values considered. 
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Table 1: Estimates for the exponents of the power functions describing the empirical dis- 
tribution tail decay for log returns over periods of different lengths 5. Here f3± are obtained 
as the slope coefficients of the fitted straight lines in Fig. [ij while /3±,{p} are Hill estimators 
based on the left-most/right-most p% of the data. 
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2.2 Rank correlations 



Now we will turn to analysing the dependence structure of our time series. The stan- 
dard approach for such a task is based on working with linear correlation coefficients. 



However, as we saw in Section 2.1, we are dealing here with heavy-tailed data for 
which using linear correlations may be unwise. For that reason, and also because our 
intention is to use for modelling purposes a transformation of a Gaussian process with 
a simple dependence structure, we will prefer to analyse rank correlation coefficients in 
this section (as they are invariant under strictly increasing transformations and thus 
may be useful for doing statistical inference for the underlying unobservable process). 

As is well known, the values of the two most popular rank correlation measures, 
Kendall's tau and Spearman's rho, are quite close to each other (which was also con- 
firmed by the data we deal with in this section). Because computing the latter requires 
less computational effort, we chose to work with Spearman's rho. 

Recall (see e.g. [18], pp.207, 229) that the theoretical value of Spearman's rank 
correlation coefficient ps for a pair of random variables U, V is defined by 

Ps (U,V)=p(Fu(U),F v (V)), (2) 

where Fjj and Fy are the (marginal) distribution functions of U and V, resp., and p is 
the linear correlation. The standard estimator of ps{U, V) from a sample (Uj, Vj), j = 
1, . . . , n, is calculated using the ranks of the variables within the respective univariate 
samples as 

*> ■= t (™ k w> - n -V) - n -V) ■ P> 

It is clear that calculating rank correlations for the whole sample of daily data over 
the eight year long period, without careful detrending and possibly some further data 
pre-processing, would hardly be meaningful. Furthermore, as we pointed out earlier, 
we would be most interested in modelling the price processes over medium-term time 
intervals (about 100 days long), over which one can expect data to follow a reasonably 
stationary process. And indeed, computing correlations for data from time intervals 
of such lengths leads to rather stable results. A typical representative is depicted in 
Fig. [3j showing Spearman's rho values for the data consisting of the pairs (Z t , Z t+ k), 
t — 0,1, ... , 200 — k, for values k = 1, . . . , 35. One can observe that the plotted line 
looks pretty much like an exponential curve, of the form e~ ck , a behaviour typical of 
autoregressive processes of the first order (or OU processes in continuous time). 

For the increments of the time series values, i.e. the log returns over time periods 
of a fixed length 5, the picture is even more interesting. Following our general ap- 
proach, here we will pool data from the whole eight year long period of observations, 
as the influence of the possible (slow) trend present in the time series of the values 
of the increments of Z t over relatively short lags would be quite small. Fig. [4] shows 
three plots of Spearman's rho's values calculated for the data consisting of the pairs 
(AgZ t , AgZ t+ k), t = 5, 5 + 1, . . . , N — k, that are plotted for k — 1, 2, ... , 100, for the 
values of 6 equal to 20, 40 and 60. 
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Figure 3: Spearman's rho for (Z t , Z t +k), k = 1,.. . , 35, from the first 200 data points. 
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Figure 4: The plots of Spearman's rhos for the pairs of log returns (A$Zt, A$Z t +k), t = 
6, 8 + 1, . . . , N — k, over time periods of 5 = 20, 40 and 60 days, plotted for k = 1, 2, . . . , 100. 

A remarkable common feature of these plots is that they decay in an almost linear 
way for lag values k from 1 to 5, attaining (small) negative values at the minimum 
points, and then start slowly growing. 



2.3 Copulas 

To further analyse the dependence structure of our data, we turn to copulas. One 
standard way of graphical representation of empirical copulas is to make scatterplots 
of the pseudo-sample obtained by transforming the components of the sample points 
using the marginal empirical distribution functions (see e.g. p. 232 in [T8]). 

We will begin with considering the original data set {Z t }. To avoid the interfer- 
ence of the possible long-term trend, first we process the original data from blocks of 
length 100. We consider 20 subsamples 

{(U ktj , V kJ ) := (Z tk+j , Z tk+j+5 ), j = 0, 1, ... , 94}, where t k = 100k, k = 0, 1, . . . , 19, 
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and for each of them produce a sub-pseudosample given by 



i = <>,..., 94, 

where rank(C4j) is the rank of Uk,j in the fcth subsample Uk,o,Uk,i, ■ ■ ■ ,^,94, and 
likewise for rank(Vfcj). The aggregate pseudosample is then obtained as the union of 
these 20 sub-pseudosamples and is depicted on the left pane in Fig. |5j 




Figure 5: The aggregate scatterplot for the pseudosamples obtained from the data (Z t , Zt+5) 
(left pane) and the scatterplot of N simulated points from the Gauss copula C^ a with p = 
0.693 obtained from the estimate of Spearman's rho. 

Next to that scatterplot representing the dependence structure for the log prices at 
lag five, we put the scatterplot of N points simulated from the Gauss copula 

$ being the standard normal distribution function and $ _1 its inverse, with the pa- 
rameter p value obtained from our estimate for Spearman's rho. Recall that, for any 
bivariate random vector (Xi^X?) following a meta-Gaussian distribution with cop- 
ula Q, one has 

6 p 

p s (X 1 ,X 2 ) = - arcsin - (5) 

7T Z 

(see e.g. Theorem 5.36 in [18]), so that we can use the estimate p = 2sm(nps/6) for p. 
In the case of the data set (Z t , Z t+5 ), this leads to p — 0.693. 

No localisation of data followed by the subsequent aggregation is needed for the 
increments of the original time series values. Figure [6] presents the scatterplot of the 
pseudosample constructed from the data set {(A 20 Zt, A 20 Zt+5), t — 20, 21, . . . , N — 5}, 
alongside with the scatterplot of a sample simulated from the Gauss copula with the 
parameter value p = 0.695 (also estimated via (JsJ) from the respective Spearman's rho). 

We can see that, in both cases, there is a reasonable agreement of the empiri- 
cal data with the fitted Gauss copulas. To confirm this observation numerically, we 
used a recently suggested method for calculating approximate P-values for testing the 
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Figure 6: The scatterplot for the pseudosample obtained from the data (A20 Zt, /S.2oZ t +§) (left 
pane) and the scatterplot of N simulated points from the Gauss copula Cp with p = 0.695 
obtained from the estimate of Spearman's rho. 

goodness-of-fit by parametric copula families (T5J HE] . As the method is valid for i.i.d. 
samples, to use it in our situation we first had to weaken the dependence between sam- 
ple points by "rarefying" the data. We applied the test to samples obtained by taking 
each 10th pair of log-returns over two time intervals of length 5=1,5, 10 and 20 days, 
one shifted by 5 days relative to the other, thus dealing with four bivariate samples 
each consisting of 206 points of the form (AgZ t , AsZ t+ z,), t = 10, 20, 30, ... , 2060. The 
method uses Monte Carlo techniques, and based on 100 simulations from the respective 
Gauss copulas for each of the 5 values considered, we obtained the following P-values: 
0.675 (5 = 1); 0.332 (5 = 5); 0.214 (5 = 10); 0.433 (5 = 20); so that the Gauss copula 
hypothesis was not to be rejected. 

In conclusion of this section, we summarise its key findings that will be referred to 
in Section [3l 

[Fi] The log returns A$Z t over periods of 5 days have distributions of which the 
"central parts" (about 95% of all data) are very well fitted by the respective 
parts of normal distributions. 

[F2] The empirical distributions of the log returns A$Z t have power tails, with the 
exponent of the power function (roughly) independent of the lag 5. 

[F3] Spearman's rho for (Z t , Z t+ k) decays as an exponential function of k. 

[F 4 ] As a function of k, Spearman's rho for (A$Z t , AgZ t+ k) behaves in the fashion 
presented in Fig. |4j 

[F 5 ] For fixed k and 5, the empirical copulas for both (Z t , Z t+k ) and (AgZ t , A$Z t+k ) 
agree with the respective Gauss copulas C^ a , with the parameter p estimated 
from the calculated values of Spearman's rho for the samples. We illustrated 
that for k = 5, 5 = 20, and gave approximate P-values for a goodness-of-fit test 
for Gauss copulas for (AgZ t , A$Z t+ k) with k = 5 and 5 = 1,5, 10 and 20, but the 
situation is similar for other values of the quantities as well. 
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3 The model and its key properties 



In this section, we present a formal description of our simple model suggested by the 
findings [Fi]-[F 5 ] and demonstrate that it does have properties consistent with these 
empirical facts. 

Let {X t } be a stationary OU process driven by the stochastic differential equa- 
tion (11]). It is well known that {X t } is a Gaussian process, with X t ~ N(0, a 2 ), where 
a 2 = r 2 /(2a). Furthermore, for < s < t, the conditional distribution of X t given 
X s = x is normal with mean and variance given by m t - s (x) and c 2 _ s , respectively, 
where we used notation 
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m u (z):=ze- au , a 2 u := T — (1 - e" 2 *"), 

so that the (linear) correlation function of the process has the form 

p(X s ,X,) = e-^- s l. (6) 

For a strictly increasing continuous mapping h : R — > R, consider the process 

Y t := h(X t ), t > 0. 

It is obvious that {Y t } is also a stationary Markov process and, moreover, that if 
h is twice continuously differentiable then {Y t } will be a diffusion process as well. 
Observe that all finite-dimensional distributions of {Y t } are meta-Gaussian: for any 
< ti < t 2 < ■ ■ ■ < t n , the vector (Y tl , . . . , Y tn ) has a Gauss copula with the correlation 
matrix (p^ = e ~ a ^ ti ~ t ^) i j <n an d identical univariate marginal d.f.'s all equal to F(y) : = 
]P(Vt < y) = ^(h^ 1 (y)/a), $ being the standard normal d.f. 

Our first objective will be to determine conditions on h under which one or both of 
the distribution tails of Y t is/are of regular variation, so that for the d.f. F of Y t one 
has 

F(y):=l-F(y) = y -^L + (y) (7) 
and/or 

F{-y)=y-f t -L-{y), (8) 

where (3± > and L± are functions slowly varying as y — > oo. 

To this end, we will introduce two conditions [A + ] and [A_] as follows: 

[A-t] For some constant b± > and differentiable function f±(x), one has 

h(x) = ± exp{b±x 2 + f±(x) + o(l)} with f' ± (x) = o(x) as x — » ±oo. (9) 
Now we can state the following key result on the tail behaviour. 

Theorem 1. The following assertions hold true. 

(i) Under condition [A + ] ; the stationary distribution FofY t has a right tail of the 
form with (5+ = l/(2b + a 2 ) = a/{b + r 2 ). 



10 



(ii) Similarly, under condition [A_] the left tail of F is of the form (J8|) with /?_ = 
l/(26_a 2 ) = a/(^r 2 ). 

(iii) Under condition [A+], i/ie ng/ii toz/ of the transition distribution function 
F s , t (y s , yt) ■= W(Yt < y t \Y s = y s ), < s < t, is regularly varying of index 

_n _ 1 a _ 

P+ ^- s '- 2b + a 2 t _ s ~ 6+r 2 (l -e-Mts)) < P+ ' 

Under condition [A_], a symmetric assertion holds for the left tail of the transition 
distribution, with the regular variation index — (3-j- s '■= — (2&-0" 2 _ s ) _1 - 

(iv) If [A + ] is met and h(—x) = o{h(x)) as x — >■ oo ; then, for any < s < t, 
the right tail of the distribution of the increment Y t — Y s in the stationary process has 
the same regularly varying asymptotic behaviour at infinity as that of Y t from part (i) . 
Similarly, if [A_] holds and h(x) = o(h(—x)) as x — > oo, then the right tail of the 
increment Y t — Y s in the stationary process has the same regularly varying asymptotic 
behaviour at infinity as the left tail of F from part (ii) . 

The proof of the theorem is given in Section [5j 

Now suppose that our h is linear in the "middle part": say, on an interval / such 
that F(X t £ I) > 1 — e/2 (with a small enough e), we have 

h(x) = a + cx, x <E I, 

for some constants a G 1R, c > 0, and then outside / the function h has "tail parts" 
satisfying conditions [A±]. Then, assuming that the log price time series {Z t } is 
represented by (the values at integer time points t of) the process {Y t = h(X t )}, we 
will have for the log returns a representation of the form 

A 5 Z t = h(X t ) - h(Xt- S ) = c(X t - X t _s) + P(&i + 0) < e, (10) 

since {^ t ,s ^ 0} C {X t 1} U {X t -& /}. Therefore the "central part" of the distribu- 
tion of A$Z t will be very close to that of c(X t — Xts), which is clearly normal. This 
shows that our model is consistent with empirical fact [Fx] . 

That [F2] is also reproduced by the model follows from Theorem [l](iv). Here a 
curious observation is in order. On the one hand, the log-log plots in Fig. [2] dis- 
play "shifted" parallel straight line segments in the curves corresponding to different 
values of 5, and such a translation indicates that the tail behaviours of the returns' 
distributions for different 5's differ by constant factors (increasing with 5). On the 
other hand, Theorem [l](iv) claims that they should have common asymptotics, which 
seems to be a contradiction. However, large sample simulations of the differences of 
transformed components of Gaussian vectors demonstrate the same translation of the 
straight line segments present in the log-log plots of the empirical distribution tails for 
increments, corresponding to different 5 values, in the "moderately large" deviations 
zone (in agreement with the empirical observations), which is then followed by a "fu- 
sion" of the curves for larger deviation values thus confirming the common asymptotics 
established in Theorem [l|iv). The proof of the theorem suggests the following expla- 
nation for this phenomenon: the established (common) asymptotics are essentially due 
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to the convergence of the conditional probability in the integral on the right-hand side 



of (21) (or, rather, the conditional probability f(Y s < z — y \ Y t = z) in ^3, see (16)) to 
one. However, in vicinity of the lower integration limit (1 + e)y in ^3 this convergence 
is slow. In combination with the regular variation of the right tail of the distribution 
of Y t , this results in a "pre-limiting" tail beahaviour differing from that of ¥(Y t > y) 
by a factor depending on 6 = t — s and y, which can be small for small 6 (hence the 
translation of the plots), but which eventually tends to one as y — > 00. 

To illustrate the above statement concerning the simulation study, Fig. [7] displays 
the log-log plots of the empirical distribution tails for three i.i.d. samples (of size 3 x 10 5 
each) of the differences h(( 2 ) — MCi); where 

h(x) =2(e {3+sisnx)x2/20 -l)signx + x, xeR, (11) 

has the form stipulated in ([9]), with b + = 0.2 and = 0.1, and (£1, £ 2 ) have standard 
normal bivariate distributions with correlations e~ a , with a = 0.5 (the left curve), 1 
and 1.5 (the right curve), respectively, playing the role of a8. 




Figure 7: The log-log plots of the empirical distribution tails for simulated differences 
MC2) — ^(Ci) f° r i-i.d. standard Gaussian random vectors ((1,(2) with different correlations 
(for detailed description, see the text around formula (|11|)). 



That property [F 3 ] holds for {Z t := h(X t )} follows from the exponential form of 
the correlation function ([6]), the invariance of Spearman's rho under strictly increasing 
transformations (so that ps(h(X s ), h(X t )) = p s (X s , X t )) and relation ^ (note that 
the function of p on the RHS of ^ is very close to the identity function on the interval 

pe[-l,l]). 

The demonstration of [F4] is somewhat less straightforward. First we state the 
following bound of which the proof is given in Section |5j As can be seen from the 
proof, the bound is rather conservative, and so the result is more qualitative than 
quantitative in nature: one can expect that, in most cases, the distance between the 
values of the Spearman's rhos will be much smaller than the bound given. 

Theorem 2. Let U , V , £ and rj be random variables given on a common probability 
space and such that U, V and the sums 
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are all continuously distributed. Then, setting P(£ 7^ 0) =: E\ and P(r/ 7^ 0) =: 62, one 
has 

\ps{Ui, V x ) - Ps {U, V)\ < 18(ei + e 2 ) + 12e l£2 . 



From that bound and ( J10[ ) it follows that 

\p s (A 5 Z t , A s Z t+k ) - p s (A s X t , A s X t+k )\ < e , 

where e = 36e + 12e 2 . Here (A 5 X t , A^Xt+fc) = (X t - X t ^ 5 , X t+k - X t+k ^ s ) has a joint 
Gaussian distribution, and so in view of the above-stated relationship ^ between ps 
and p in the Gaussian case, it remains to evaluate the linear correlation 

i\ v a v \ ^( X t - X t - S )(X t+k - X t+k - S ) 

P[^S^t, ^S^t+k) = 



Xar(X t -X t _ s ) 

where we used the stationarity and zero means of {X}. Expanding the numerator and 
employing ([6]), we obtain, for k > and 5 > 0, 

E(X t — X t -5){X t+k — X t+k s) 

= KX t X t+k — EX t X t+k -s — E^i-^i+i; + EX t _sX t + k s 
= a 2 (2p(X t , X t+k ) - p{X t , X t+k _ 5 ) - p(X t , X t+k+5 )) 
= o- 2 (2e- ak - e~ alk - 51 - e- a{k+5) ). 

As a special case of this where k = 0, we have Var(X t — X t -s) = 2a 2 (l — e~ aS ), so that 

2e~ ak — e _a ' fc_5 ' — e~ a ( k+s ^ 



p(A 5 X t , A S X, 



t+kj 



2(1 -e~ aS ) 

The plots of this function of k for three different values of S are presented in Fig. [8j 




Figure 8: The plots of the correlation functions p(A$Xt, AsX t -t- k ) in the case a = 0.05, for 
5 = 20, 40 and 60 days, plotted for k < 100 (cf. Fig. Q. 

They demonstrate virtually the same behaviour as the empirical Spearman's rhos shown 
in Fig. [1J thus confirming that our model has property [F 4 ] . 



13 



Concerning [F 5 ], we immediately see that, for our model Y t = h(X t ), the vector 
(Yt,Y t+ k) clearly has a Gauss copula. As it was the case when dealing with [F 4 ], 
showing that the copula for the vector of increments (AgYt, AsY t+ k) will be close to a 



Gauss one, provided that the value of e in (10) is small, is a more sophisticated task. 
However, that can be easily done adapting the proof of Theorem [2] and its extension to 
the case of bivariate distribution functions. This is a straightforward technical exercise 
that we leave to the interested reader. 

Remark 1. A possible relatively simple modelling approach alternative to ours would 
be to use, instead of (jl]), a stochastic driver of the form 

dX t = -aX t dt + dL t , (12) 

where {L t } is a suitably chosen Levy process. In particular, Theorem 3.2 in [9] estab- 
lishes the following result: Let /j el and u, A > be arbitrary, 

If 00 r x ~\ 

Kv(x) = -J y' / - 1 exp{--(y + y- 1 )jdy 

1 f°° 

= — I exp \yu — x coshw} du, x > 0, 

* J -co 

be the modified Bessel function of the third kind, and {L t }t>o be a Levy process with 
cumulant function 

lnEe» Ll = iun - \\ u \ ^fi=£^ , ueR\ {0}. 

K u/2 {X\u\) 



Then the stochastic differential OU type equation (12) has a weak strictly stationary 
solution satisfying 



X t = e~ at X + e~ at I e as dL s 
Jo 

and such that X t has marginal t- distribution T{y, A,/i) with density 

c(v, A) 



x G 



where c(v, A) = (AB(z//2, 1/2)) -1 , B(-, •) being the beta function. Moreover, if v > 1 
then E,X t = fi, and if v > 2 then corr (X t+T , X t ) = e~ a ' T ' . 

As the resulting process {X t } will already have power function decay of the distri- 
bution tails, it might be possible to use that process for modelling the financial data of 
interest, without transforming it (although some transformation might still be needed, 
due to the symmetry of t distributions). The advantage of using a transformed OU 
model of the type dealt with in the present paper is that it is more flexible as its 
temporal dependence structure can be decoupled from its distribution tail properties. 
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4 Further comments on model fitting 



First we observe that, without loss of generality, one can always put r = 1 in 
indeed, we clearly have Y t = h(X t ) with h(x) := h{rx) and X t := T~ x X t satisfying 

dX t = -aX t dt + dW t , t > 0. 

So in the "driver" component Q of the model, one only needs to estimate a, and it is 
this task that will be discussed in the present section. 

Concerning the estimation of the function h, we will only point out that, in view of 
its specific form dictated by Theorem [TJ a parametric approach would be most natural. 
One possible simple parametric class of functions consists of candidates of the following 
spline type: 

h(x) = —aie b ~ x2 l(x < X\) + (a 2 + a 3 x)l(xx < x < x 2 ) + a^^Xix > x 2 ), 

where l(-) denotes the indicator function, the bounds Xi can be estimated as the end 
points of the "central part" of the distribution well approximated by a normal one, a 2 
and (23 > can be estimated from that "central part" (i.e. from a truncated sample), 
parameters b± > can be estimated using estimates for the tail regular variation indices 
(e.g. Hill estimators) and the assertions of Theorem [l] (note that, to use them, one has 
to estimate a first), and a\ and 04 can be chosen to make the function h continuous. 
Alternative parametric families may be more adequate. 

Suppose that the process {Y t ,t G [0, T]} is observed at times t k = ke, k = 0, . . . , n, 
where n + 1 is the total number of observations. Without loss of generality, we can 
always assume that e = 1. 



4.1 An estimator based on rank correlations 

Since Spearman's rho is invariant under strictly increasing transformations of data and 
(Xj,Xk) is Gaussian, we have from ^ that 

p s {Yj, Y k ) = p s (Xj,X k ) = - arcsm = - arcsm , 

so that 

a = -^ ln ( 2sm & s(y ' ,n) ))' 

Hence one can use the "plug-in" estimators 

a k := -iln(2sinQp s (r ,n))) , k = l,2,..., 

where ps(Yo,Yk) is calculated according to (|3]), for the sample of pairs (Y t ,Y t+ k). As 
discussed before, such an estimation would require blocking of data, to reduce the effect 
of the slow trend component in the time series. Alternatively, one can use Spearman's 
rho for the increments of Y t (cf. our discussion following Theorem [2]). The values of a k 
for different fc's obtained from 50 day log-returns for BHP data are shown in Table [2] 
Note the stability of the estimates, which confirms the appropriateness of the model 
(according to which the quantities are estimators of the same value). 
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k 


5 


10 


15 


20 


25 


30 




0.032 


0.036 


0.036 


0.037 


0.036 


0.038 



Table 2: The values of from 50 day BHP stock log returns, for different values of k. 



4.2 An alternative estimator 

An alternative (but closely related) estimator for a is based on the following simple 
observation concerning the background OU process, which we will formulate in the 
form of a theorem (for its proof, see Section [5]). Set 



g{x) := - In [sin (2tt(x - 1/4))] , x G (1/4, 1/2). 

Theorem 3. Let 



I n—k 

p2:=^E 1 ( A i >0 ' A >+* > °) 



be an empirical estimate of p^ := P(X > 0,X k > 0). Then 

& k := g -^, k = l,2,..., (13) 
are strongly consistent and asymptotically normal (as n — > oo) estimators of a with 



Var(a 



k) ^ 



vl(9'(Pk)) 



2 



nk z 
where 



v\ : = p fc (l - Pk ) + 2J2 m*k+i > 0, X, > 0, X fc > 0, X > 0) - pi) . (14) 



z=i 



Of course, in the context of our problem {X t } is unobservable, so we need to replace 
p° k with an estimator for p k from {Y t }. Denoting by m the median of F, we clearly 
have {Xj > 0} = {Yj > m} and so 



n— fc 



We still do not know m, but one can use instead the empirical median rh n (Y) of the 
sample Y ,Yi, . . . ,Y n . It is not hard to verify that m n {Y) m as n — > oo, which 
implies that 

. n—k 

Pk ■= 7 ^2 > "^(F), ^j+fc > ™>n{Y)) Pfe- 

i=o 

Therefore 

- ._ ?M k _ 

a k .— , At — i, z, ... , 
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k 


5 


10 


15 


20 


25 


30 




0.027 


0.035 


0.035 


0.032 


0.031 


0.031 



Table 3: The values of a' k from 50 day BHP stock log returns, for different values of k. 



will also be strongly consistent estimators of a. 

The values of a' k for different k's obtained from 50 day log-returns for BHP data 
are shown in Table HI 

Simulations showed that the root mean square errors of the estimators pk are smaller 
than those for p° k , so that from the viewpoint of the quadratic error minimisation 
it is better, in a sense, not to know the true value of the median of the stationary 
distribution F. Figure [9] displays the values of the ratio of the standard error of p k to that 
of pk, for k = 20, calculated from a sample of 10 4 independent discretised trajectories 
of the OU process on the time interval [0,T], with T = 10 4 and discretisation step 
e = l. 




2.5 ' ' ' ' ' ' ' ' ' 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 9: The ratio of the root mean square errors of P20 ano ^ P20 as a function of a 
calculated using crude Monte-Carlo simulations of the respective OU process. 

It can be seen from the plot that, as a increases, the value of the ratio appears 
to decrease to some limiting value. That limiting value should correspond to the case 
of independent observations Yj, j = 0,1,2, ... ,n (or Xj), in which the ratio can be 
calculated, as shows Theorem [4] below (its proof is given in Section [5]). Unfortunately, 
we were not able to quantify the effect in the general case, for finite values of a, but 
a plausible "common sense explanation" of it remains the same: the "quality" of the 
estimators pi and pk is basically determined by how often the difference Xj — a changes 
its sign, for a = and a = rh n (X), respectively, and in the latter case this occurs more 
often. 

For notational convenience, in Theorem [4] we will consider samples starting from 
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Xi (rather than X ), so that here we will have 

- n—k 

p2:=^rE 1 (- X 'i >0 '- X i+* >0 )' 

3=1 

Pfc := > m n (X),X i+fc > m n (X)), 

n — k 

3=1 

m n (X) being the sample median for X\, X%, . . . , X n . 

Theorem 4. Let Xi,X%, . . . ,X n be i.i.d. standard normal random variables and p k '■ = 
P(Xj > 0,Xj + k > 0), k = 1, 2, . . . Then, for any fixed k, p Q k is an unbiased estimator 
°f Pk = 1/4, pk has a bias of 0{n~ l ) as n — >■ oo ; and 

Var (p°) = -|- + 0(n- 2 ), Var (p fc ) = E(p k - p k f + 0(n- 2 ) = -L + 0(n- 2 ). 



4.3 An estimator based on band crossing times 

In the two previous subsections we considered estimation of the drift coefficient a of the 
unobservable OU process {X t } from its transformed version {Y t = h(X t )} which was 
assumed to be observable. It is interesting to note that one can estimate a even if {Y t } 
is only partially observable itself — namely, if we only know the number of times the 
process' trajectory crossed the band between two given levels A and B during a given 
time interval, provided that we also know the values of the d.f. F(y) of the stationary 
distribution of Y t at the points y = A, B. 

For a, b E E denote by T a b the time it takes the OU process {X t } driven by ([T]) 
with r = 1 to transit from point b to point a: slightly abusing notation, one can write 

T afi := M{t >0:X t = a\X = b}. 

The density of T a ^ has a simple closed-form representation when a = [2E]; in the 
general case, this is not so (one needs to invert a Laplace transform expressed in terms 
of parabolic cylinder functions, see e.g. [21] )• However, for our purposes it will be 
enough to know the means of T a ^ only, and these are available, see [24J . 

By the excursion time T( a ),b to a from b we will call the time it takes the OU process 
starting from b to reach the point a and then return back to b. Clearly, 

T( a ),b — T a ,b + Tb, a , 

where the random variables on the right-hand side are assumed to be independent. 

Set a := /i _1 (v4), b := h^i^B). Since X t ~ iV(0, l/(2a)) in the stationary regime, 
one has §{^/2aa) = F(X t < a) = T>(h{X t ) < A) — F(A) and likewise §(V2ab) = 
F(B). Therefore 

a _<$>- 1 (F(A)) a* b _^~ 1 (F(B)) _ b* 
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The standard scaling argument shows that the process X* := a 1/2 X t/a satisfies Q 
with a = 1, r = 1 (and another Brownian motion process, namely W£ = ot l l 2 Wt/ai 
instead of W-t). Therefore one has 

where the asterisk indicates that the excursion time on the right hand side is for the 
process {Xf}, the distribution of T? a *-> b * being independent of a and having a mean 
iTL a *,b* which can be computed from F(A) and F(B) using the results of |24j . 

Now the strong Markov property implies that the subsequent excursions from b to 
a in {X t } (or, equivalently, from B to A in {Y t }) form an i.i.d. sequence. Therefore, 
if N^b{T) denotes the number of times {Y t } crosses the band between the levels A 
and B in a given direction (say, from B to A, each such crossing corresponding to an 
excursion from B to A) during the time interval [0, T], the integral renewal theorem 
asserts that, as T — > oo, 

Na,b( t ) a.s. 1 a 
T ET (a))6 m a .,6*' 

which leads to another consistent estimator of a given by a" := m a *^*N^B{T)T~ l . In 
conclusion we note that one can evaluate the mean quadratic error of the estimator 
a" from the central limit theorem for renewal processes (the variance of T a ^ is also 
available from [21]). 



5 Proofs 

Proof of Theorem^ (i), (ii) As part (ii) is a "mirror reflection" of part (i), we only 
need to prove the latter (in fact, part (ii) was included in the theorem as a separate 
assertion only for convenience of referencing to it in part (iv)). 

Set x := h^ 1 {y). Since X t ~ iV(0,cr 2 ), we can use Mills' ratio for the normal tail 
and the relation b + /3 + = l/(2o~ 2 ) to write that, as y — > oo, 

¥(Y t >y) = F(X t >x) = -^— e - x2/{2a2 \l + o(l)) 

XV27T 



" f e b + x2 ) P+ (l + o(l)) = ^=e^^ + °( 1 \h(x))-P+ 
x\/2vr ^ ' xy2TC 

a e 0+f + (x)-inx+o(l) y -0+ _ ( 15 ) 



2vr 

It remains to show that g(y) := e^ x \ where fo(x) = (3 + f + (x) — Inx and x = h^ 1 (y), 
is a slowly varying function of y. To this end, fixing an arbitrary A > 1 and setting 
s := /i _1 (A?/) — x, we can write 

^TT = ex vM x + s )~ fo(x)}. 

g(y) 

Since clearly fo(x) = o(x) as x — )• oo, it now suffices to show that s = 0{l/x) . 
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From the definition of s, the mean value theorem and the assumption f' + (x) = o(x), 
we have 

A = ~h7~V^ = exp{26+xs + b + s 2 + f+(x + s) - f+(x) + o(l)} 
h(x) 

= exp{26+xs(l + o(l)) + bs 2 + o(l)}, 

which immediately implies that s = 0(l/x) and hence completes the proof of part (i). 

(iii) We will only consider the case where [A + ] holds true. The case of [A_] is dealt 
with in the same way. 

First observe that, given Y s = y s (or, equivalently, X t = x s := h~ 1 (y s )), we can 
write the increment of the process {Yt} on [s,t] as 

Y t -y s = h(X t ) - h(x s ) = hi(X t - m t - s (x s )) - h(x s ), 

where, setting for brevity m := mt- s (x s ), the function 

hx(x) = h(x + m) = exp{bx 2 + 2bmx + m 2 + f(x + m) + o(l)} =: e k2+/l(l)+o(1) 

is again of the form (J9j), with fi(x) := 2bmx + m 2 + f(x + rn) clearly satisfying the 
relation f[(x) = o{x). Since, as we have already recalled, the conditional distribution 
of X t — mt- s (x s ) given X t = x s is N(0, u 2 _ s ), the proof of (iii) is completed by the same 
argument as used to demonstrate the first part of the theorem. 

(iv) We will prove the first part of the assertion, assuming that [A + ] is met and 
h(—x) = o(h(x)) as x — >• oo. The second half will then immediately follow, in view of 

the observations that Y t — Y s = —(Y s — Y t ) and (Y s ,Y t ) — (Y t ,Y s ), the latter being a 

consequence of the obvious relation (X s , X t ) — (X t , X s ) (for the stationary OU process, 
of course). 

Letting y — > oo and e = e(y) I slowly enough (so that, in particular, ey — > oo; 
further conditions on the behaviour of e will appear later), we can write 



¥(Y t -Y s >y) = j¥(Y s <z- y, Y t e dz) 



(l-e)j/ f(l+e)v poo 

+ =:h + I 2 + h- (16) 

(l-e)y J(l+e)y 



First consider 

.(l-e)l/ r(l-e)y 



F(Y s <z-y,Y t edz)< / P(Ys < —ey, Y t G dz) 

■oo J — oo 

< -ey, Y t <(l- e)y) < F(Y S < -ey) = F(-h(X s ) > ey). (17) 



In the case where inf x h(x) > — oo one clearly has I\ — for all large enough y. In 
the alternative case, the function v(w) := —h~ 1 (—w) will be defined for all w G R 
and have the property v(w) — > oo as w — > oo. Moreover, due to the assumption that 
h(—x) = o(h(x)) as x — > oo, there exists a number xq G IR and a continuous function 
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i](x) > vanishing at infinity such that h(—x) = —i](x)h(x) for all x > Xq] note that 
77* (x) := sup u>a . rj(u) j. as x — > 00. From that and the symmetry of the distribution 
of X s we see that, for all large enough y, the last probability in (17) is equal to 

F{-h{-X s ) > ey,X s > v{ey)) = F(r){X s )h(X s ) > ey, X s > v(ey)) 

<F( V *{v{ey))h{X s )>ey), 

and hence 



h<F(h(X s )> y 

Now since we can assume that e — > so slowly that e/r]*(v(ey)) — > 00 (as such a choice 
is clearly possible), it follows from part (i) that Ii = o(F(y)). 

Further, it is obvious from the Uniform Convergence Theorem for slowly varying 
functions (see e.g. Theorem 1.1.2 in [8]) that 

((1 ± e)y)-^L + {{\ ± e)y) = (1 + o(l))y-^L + (y), (18) 

and so from (i) one immediately obtains the bound 

h < F(Y t e ((1 - e)y, (1 + e)y\) = F((l - e)y) - F((l + e)y) 

= (1 + o(l))y-^L+(y) - (1 + o{l))y-^L + {y) = o(F(y)). (19) 

It remains to evaluate J 3 . To this end, we note that, on the one hand, in view of 
(18) and the result of part (i), one has 

POO 

h < / HYt e dz) = F((l + e)y) = (1 + o{l))y^+L + {y), (20) 

J(l+e)y 

while on the other hand, as 1 — y/z > e/(l + e) for z > (1 + s)y, 

POO 

h> / F{Y s /z<e/{l + e)\Y t = z)F(Y t edz). (21) 

J(l+e)y 

Setting u := h~ 1 {z) and recalling that the conditional distribution of X s given Y t = z 
(or, equivalently, X t = u) coincides with the law of W :— we - "^ - ^ + a t - s Z, Z ~ 
N(0, 1), we see that the conditional probability in the last integral is equal to 



P b + (W 2 - it ) + f + (W) - f+(u) + 9{W) - 6{u) < In 



1 + e 

P( b + {2ue' a{t - s) a t ^ s Z + a 2 t _ s Z 2 ) + r(W, u) < b + u 2 (l - e - 2a{t - s) ) + In 



l + e 



where 6{u) = o(l) as u -> 00, r(W,u) := f + (W) - f + (u) + 6{W) - 6{u). Recalling 
that f' + {x) = o(x) at infinity, it is obvious that the last probability tends to one 
uniformly in z > (1 + e)y as y — > 00, provided that we choose e — > slowly enough 
(ensuring that a/| hxe\ h~ l {y)). This, together with (20) and (21), shows that 
^3 = (1 + o(l))?/~ /3+ L + (y) and so, in view of the earlier bounds for Ii and I2, completes 
the proof of the theorem. □ 
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Proof of Theorem^ Since, for a continuous random variable X, Fx{X) is uniformly 
distributed on (0, 1), we see from @ that 

/><?([/!, Fx) = 12(1^(^)^(^0 -0.25). (22) 

Now, using the standard argument we see that, for any 

F Ux {x) > F(U + f < x; £ = 0) > F(U < x) - P(f ^ 0) > F v (x) - e x , 
F Vl (x) < F(U + £ < x; £ = 0) + P(£ ^ 0) < F(U < x) + s x = F v (x) + e u 

and similarly for Fy^x), so that 

sup \F Ul (x) - Fu(x)\ < ei, sup \F Vl (x) - F v (x)\ < e 2 - 

X X 

Therefore, using 1(A) for the indicator function of the event A, we have 

Fu^UJFv^) < Fu^UJFv^VMZ = rj = 0) + 1(^ 0) + 1( V ^ 0) 
= Fu^FyMUZ = v = 0) + l(^0) + 1(77 0) 
< (i^f/) + e^FviV) + e 2 ) + l(f ^ 0) + 1(»7 ^ 0), 

so that 

EF ai (C/0^yi(F) < E(Fy(C0 + e x ){F v {V) + e 2 ) + £l + e 2 
< EFu(U)F v (V) + -(ei + e 2 ) + e 1 e 2 , 

where we used the simple fact that EF[/(C7) = EF V (V) = 1/2. This, together with a 
similar lower bound (that can be obtained in exactly the same way) and relation (22), 
completes the proof of the theorem. □ 



Proof of Theorem [3| Estimator ( 13 ) is obtained as a method-of-substitution ( "plug- 
in") estimator from the relation a = g{pk)/k which is established by inverting 

Pk = \ + 1 L arCSin ( e ~ Q ") G (1/4 ' 1/2) ' k ~ (23) 

the latter representation following e.g. from Theorem 5.35 in [15] . Since p k pt as 
n — > oo due to the ergodicity of the OU process (see e.g. Theorem 5.6 in [13j), the 
consistency of a° k follows immediately from the continuity of g (see e.g. Theorem 1 in 
Section 1.5 of [7j). 

The OU process is a-mixing with exponential rate [221 121] • Therefore {Zj : = 
l(Xj > 0,Xj + k > 0)} is also a-mixing with the same rate, and so we can apply the 
central limit theorem for a-mixing sequences (see e.g. Theorem 27.4 in [5]) to claim 
that the distribution of n l / 2 {p Q k — p k ) converges weakly to N(Q,v k 1 ) with 



v k 



E(Z - p k f + 2 HZ ~ Pk)(Zi - v k ) 
i=i 

oo 

V k (l-Vk)+2j2(^Z Z l -vl). 



i=i 
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Since EZ Z, = P(X > 0,X k > 0,X t > 0,X l+k > 0), this yields representation flli] ). 

Now as g is differentiable, g' ^ 0, the standard application of the delta-method (see 
e.g. Theorem 3 in Section 1.5 of [7j) establishes that 

Vn(9(Pk) - 9(Pk))/k = \/n(a° k - a)/k 

is asymptotically normal with zero mean and variance vl(g'(p k )) 2 k~ 2 . The theorem is 
proved. □ 

Proof of Theorem^ That p k is an unbiased estimator of p k is obvious. To compute 
the variance of p k , we set n k := n — k and lj := > 0) for brevity and note that 

t n k n k~^ n k 

i=l i=l j=«+l 

i /n k -k n k n k -l n k 

4 EE-+ E £■■■)• (*> 



A: 



i=l j=i+l i=n k — fe+1 J=i+1 



For each i in the first double sum in the last line the value j = i + k will be present 
in the inner sum, and it is clear that the corresponding term will be equal to 

Elili + felj + A;li + 2fc = 1/8, 

while all the remaining terms will be equal to 1/16. In the second double sum, one 
cannot have j = i + k, and so all the terms in it will be 1/16. Therefore 

\ i=l 7 i=n k -k+\ 

1 _2_ / y> n fc - i n k - 
~ An k n 2 \ 16 16 

~ 16 + 16n^ ~ 8nf' ^ ^ 

so that Var (p°) = E(p°) 2 - V 16 = 5/(16n fc ) - k/(8n 2 k ). 

Now we turn to p^. Recall that Xi jTl > • • • > X nj „ denote the order statistics for our 
sample and let X := (X 1>n , . . . , X n ^ n ), ) el n Assume for definiteness 

that n is odd and set v := n/2 + 1/2, so that m n (X) = X VjTl . Then by symmetry and 
the total probability formula one has 



E P* = — fl F ( X J > rh n {X),X j+k > m n {X)) 

71 t. ■ * 



n k *: 



F(X l >m n {X),X 2 >m n {X)) 
/ P(Xi > x v , X 2 > x v | X = x)¥(X E dx), 
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where, again by virtue of symmetry, for x with x\ > x 2 > ■ ■ ■ > x r , 



P 2 : = P(X X > x u , X 2 > x v I X = x) 

= P ({Xi, X 2 j C {X l n , X 2i n, . . . , X u _i )n }) 

__ (Y) _ U ~ 3 







4n 



so that 



Ep fc = P 2 = l/4-3/(4n). 



(26) 



(27) 



To find the second moment of pk, we use a representation similar to (24), the total 
probability formula and a calculation similar to (25) to conclude that (keeping k fixed) 



9 1 2 

Ep 2 , = —P 2 + -3 
n fc n£ 



(tij. - z)P 4 + (n fc - fe)(P 3 - P 4 ) 



i=l 



P 4 + n- x (P 2 + 2P 3 - 3P 4 ) + 0(n 



(28) 



where, similarly to (26), 



P 3 := P(Xi > x u , X 2 > x u ,X 3 > x u \ X = x) 



and 



Cl 1 ) _ (n-3)(n-5) 



(3) 



P 4 := P(Xi > x l/ ,X 2 > x u ,X 3 > x u ,X^ > x u I X = x) 



Now substituting the values of P r into (28) we obtain that 
1 5 



16 16n 



0(n" 



Var(p fc )=Ep 2 -P 2 



3n(n - 2) 



(V) = (n-5)(n-7) 
(2) 16n(n - 2) 



+ 0(n" 2 ), 



16n 



and the quadratic error E(j3fc — p^) 2 will have the latter form, too, since Kpk — Pk = 
0(pr l ). 

If n is even, setting f := n/2 + 1 leads to the same expressions for P r in terms of v 
and n, resulting in the values 



P, 



n-2 



1 1 



4(n-l) 4 4n 



+ 0(n- 2 ), 



P 



n - 4 
ifn- 1^ 



(n-4)(n-6) 
16(n-l)(n-3)' 



so that 



+ 0(n _ 



16 16n 

as before, which completes the proof. 



Var(p fc )=Ep 2 -P 2 2 



16n 



+ 0(tT 



□ 
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